classdef EddyStatisticClass < handle
    %EDDYCLASS 此处显示有关此类的摘要
    %   此类用于涡旋的统计工作
    
    properties
        AEEDC
        CEEDC
        extent
        interval
        mask
    end
    
    methods
        function obj = EddyStatisticClass(AEEDC, CEEDC, varargin)
            % AEEDC.eddyFile2Dict();
            % CEEDC.eddyFile2Dict();
            % AEEDC.Dict2traj();
            % CEEDC.Dict2traj();
            obj.AEEDC = AEEDC;
            obj.CEEDC = CEEDC;
            idx_extent = find(strcmp(varargin,'extent'));
            if ~isempty(idx_extent)
                obj.extent = varargin{idx_extent + 1};
            end

            idx_interval = find(strcmp(varargin,'interval'));
            if ~isempty(idx_interval)
                obj.interval = varargin{idx_interval + 1};
            end     
        end

        
        function extentSet(obj, extent)
            obj.extent = extent;
        end
        
        function intervalSet(obj, interval)
            obj.interval = interval;
        end

        function EddyLimit(obj,varargin)
            obj.AEEDC = EDCLimit(obj.AEEDC,varargin);
            obj.CEEDC = EDCLimit(obj.CEEDC,varargin);   
            function [EDC] = EDCLimit(EDC,varargin)
                varargin = varargin{1};
                idx_lifetime = find(strcmp(varargin,'lifetime'), 1);
                if ~isempty(idx_lifetime)
                    lifetime = varargin{idx_lifetime + 1};
                    EDC.traj = EDC.traj([EDC.traj.lifetime] > lifetime);
                end
                idx_amplitude = find(strcmp(varargin,'amplitude'), 1);
                if ~isempty(idx_amplitude)
                    amplitude_limit = varargin{idx_amplitude + 1};
                    select_idx = [];
                    for i = 1:length(EDC.traj)
                        mean_amplitude = mean(EDC.traj(i).amplitude);
                        if mean_amplitude >= amplitude_limit
                            select_idx = [select_idx i];
                        end
                    end
                    EDC.traj = EDC.traj(select_idx);
                end
            end
        end

        % 平均属性统计 包括涡的振幅 生命周期 半径等
        function  [result] = averAttributes_measure(obj, varargin)
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj([], obj.AEEDC.traj, -1);
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj(eddyArray, obj.CEEDC.traj, 1);

            idx_type = find(strcmp(varargin,'type'), 1);
            if ~isempty(idx_type)
                type = varargin{idx_type + 1};
                eddyArray = eddyArray([eddyArray.type] == type);
            end
            idx_amplitude = find(strcmp(varargin,'amplitude'), 1);
            if ~isempty(idx_amplitude)
                result.mean_amplitude = mean([eddyArray.amplitude]);
            end
            idx_radius = find(strcmp(varargin,'radius'), 1);
            if ~isempty(idx_radius)
                result.mean_radius = mean([eddyArray.radius]);
            end

            idx_lifetime = find(strcmp(varargin,'lifetime'), 1);
            if ~isempty(idx_lifetime)
                if isempty(idx_type)
                     CElifetime = [obj.CEEDC.traj.lifetime];
                     AElifetime = [obj.AEEDC.traj.lifetime];
                     result.mean_lifetime = mean([CElifetime AElifetime]);
                elseif type == -1
                    result.mean_lifetime = mean([obj.AEEDC.traj.lifetime]);
                elseif type == 1
                    result.mean_lifetime = mean([obj.CEEDC.traj.lifetime]);
                end
            end
        end

        % 区域数量统计 %数量 振幅 半径
        function [lon, lat, result] = areastatistic(obj, varargin)
            if isempty(obj.mask) 
                [obj.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval, varargin);    
            end
            if obj.mask.mask_lon(2) - obj.mask.mask_lon(1) ~= obj.interval
                [obj.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval, varargin);    
            end
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj([], obj.AEEDC.traj, -1);
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj(eddyArray, obj.CEEDC.traj, 1); 

            idx_type = find(strcmp(varargin,'type'), 1);
            if ~isempty(idx_type)
                type = varargin{idx_type + 1};
                eddyArray = eddyArray([eddyArray.type] == type);
            end

            [lon, lat, result] = number_measure(obj, eddyArray);

            function [lon, lat, result] = number_measure(obj, eddyArray)
                latitude = [eddyArray.latitude];
                longitude = [eddyArray.longitude];
                lon = obj.mask.mask_lon;
                lat = obj.mask.mask_lat;
                for i = 1:length(lon)
                    for j = 1:length(lat)
                        if obj.mask.mask(i,j) == 0
                             result.number(i,j) = nan;
                             result.amplitude(i,j) = nan;
                             result.radius(i,j) = nan;
                             result.speed(i,j) = nan;
                            continue;
                        end
                        inside_extent = [lat(j) - 1/2*obj.interval, lat(j) + 1/2*obj.interval;lon(i) - 1/2*obj.interval, lon(i) + 1/2*obj.interval]; 
                        [xbox, ybox] = EddyStatisticClass.area_box_set(inside_extent);
                        in = inpolygon(longitude, latitude, xbox, ybox);
                        eddys = eddyArray(find(in==1)); 
                        result.number(i,j) = sum(in);
                        result.amplitude(i,j) = mean([eddys.amplitude]);
                        result.radius(i,j) = mean([eddys.radius]);
                        result.speed(i,j) = mean([eddys.speed]);
                    end
                end
            end
        end
        
        % 区域涡旋极性计算 
        function [lon, lat, polar] = eddypolarstatistic(obj, varargin)
            if isempty(obj.mask) 
                [obj.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval);    
            end
            if obj.mask.mask_lon(2) - obj.mask.mask_lon(1) ~= obj.interval
                [obj.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval);    
            end
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj([], obj.AEEDC.traj, -1);
            eddyArray = EddyDatasetClass.eddyArrayAddByTraj(eddyArray, obj.CEEDC.traj, 1);
            latitude = [eddyArray.latitude];
            longitude = [eddyArray.longitude];
            lon = obj.mask.mask_lon;
            lat = obj.mask.mask_lat;
            polar = nan(length(lon),length(lat));
            for i = 1:length(lon)
                for j = 1:length(lat)
                    if obj.mask.mask(i,j) == 0
                        polar(i,j) = nan;
                        continue;
                    end
                    inside_extent = [lat(j) - 1/2*obj.interval, lat(j) + 1/2*obj.interval;lon(i) - 1/2*obj.interval, lon(i) + 1/2*obj.interval]; 
                    [xbox, ybox] = EddyStatisticClass.area_box_set(inside_extent);
                    in = inpolygon(longitude, latitude, xbox, ybox);
                    eddys = eddyArray(in==1);
                    AEeddys = eddys([eddys.type] == -1);
                    CEeddys = eddys([eddys.type] == 1);
                    numAE = length(AEeddys);
                    numCE = length(CEeddys);
                    polar(i,j) = (numAE - numCE)/(numAE + numCE);
                end
            end
        end

        % 区域涡旋传播速度
        function [lon, lat, velocity] = eddyvelocitystatistic(obj, varargin)
            if isempty(obj.mask) 
                [obj.mask.mask_lon, obj.mask.mask_lat, obj.mask.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval);    
            end
            if obj.mask.mask_lon(2) - obj.mask.mask_lon(1) ~= obj.interval
                [obj.mask.mask_lon, obj.mask.mask_lat, obj.mask.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval);    
            end

            idx_type = find(strcmp(varargin,'type'), 1);
            if ~isempty(idx_type)
                type = varargin{idx_type + 1};
                if type == -1
                    eddyArray = eddyArrayAddVelocity(obj.AEEDC.traj, -1);
                elseif type == 1
                    eddyArray = eddyArrayAddVelocity(obj.CEEDC.traj, 1);
                end
            end

            [lon,lat, velocity.u, velocity.v] = velocity_measure(obj, eddyArray);

            % 通过已经计算好的带有速度的eddyArray计算区域涡旋传播速度
            function [lon,lat, u, v] = velocity_measure(obj, eddyArray)
                latitude = [eddyArray.latitude];
                longitude = [eddyArray.longitude];
                lon = obj.mask.mask_lon;
                lat = obj.mask.mask_lat;
                u = nan(length(lon),length(lat));
                v = nan(length(lon),length(lat));
                for i = 1:length(lon)
                    for j = 1:length(lat)
                        if obj.mask.mask(i,j) == 0
                            u(i,j) = nan;
                            v(i,j) = nan;
                            continue;
                        end
                        inside_extent = [lat(j) - 1/2*obj.interval, lat(j) + 1/2*obj.interval;lon(i) - 1/2*obj.interval, lon(i) + 1/2*obj.interval]; 
                        [xbox, ybox] = EddyStatisticClass.area_box_set(inside_extent);
                        in = inpolygon(longitude, latitude, xbox, ybox);
                        eddys = eddyArray(in==1);
                        u(i,j) = mean([eddys.u]);
                        v(i,j) = mean([eddys.v]);
                    end
                end
            end

            %计算每个轨迹的各个涡旋的传播速度并将其速度属性加入到eddyArray中
            function eddyArray = eddyArrayAddVelocity(traj, type)
                eddyArray = [];
                for trajI = 1:length(traj)
                    trajItem = traj(trajI);
                    eddys = EddyDatasetClass.eddyArrayAddByTraj([], trajItem, type);
                    eddys = eddysAddVelocity(eddys);
                    eddyArray = [eddyArray eddys];
                end
                function [eddys] = eddysAddVelocity(eddys)
                    n = length(eddys);
                    for eddysI = 1:n
                        if eddysI~=1&&eddysI~=n
                            eddys(eddysI).u = (eddys(eddysI + 1).longitude - eddys(eddysI - 1).longitude)/(2*24*3600) * 111000;
                            eddys(eddysI).v = (eddys(eddysI + 1).latitude - eddys(eddysI - 1).latitude)/(2*24*3600) * 111000;
                        end

                    end
                    eddys(1).u = eddys(2).u;
                    eddys(1).v = eddys(2).v;
                    eddys(n).u = eddys(n-1).u;
                    eddys(n).v = eddys(n-1).v;
                end
            end
        end

        
        % 计算涡旋的空间生消分布特征
        function [lon,lat,Generationstatistic,ExtinctionstatisticResult] = eddyGenerationAndExtinctionstatistic(obj, varargin)
            if isempty(obj.mask) 
                [obj.mask.mask_lon, obj.mask.mask_lat, obj.mask.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval);    
            end
            if obj.mask.mask_lon(2) - obj.mask.mask_lon(1) ~= obj.interval
                [obj.mask.mask_lon, obj.mask.mask_lat, obj.mask.mask] = EddyThermalTransortClass.createMask(obj.extent, obj.interval);    
            end


            idx_type = find(strcmp(varargin,'type'), 1);
            if ~isempty(idx_type)
                type = varargin{idx_type + 1};
                if type == -1
                    [eddyGenerationArray, eddyExtinctionArray] = GenerationAndExtinctionEddyArray(obj.AEEDC.traj, -1);
                elseif type == 1
                    [eddyGenerationArray, eddyExtinctionArray] = GenerationAndExtinctionEddyArray(obj.CEEDC.traj, 1);
                end
            end
            [~, ~, Generationstatistic] = GenerationAndExtinctionstatistic(obj, eddyGenerationArray);
            [lon, lat, ExtinctionstatisticResult] = GenerationAndExtinctionstatistic(obj, eddyExtinctionArray);

            function [lon, lat, statisticResult] = GenerationAndExtinctionstatistic(obj, eddyArray)
                latitude = [eddyArray.latitude];
                longitude = [eddyArray.longitude];
                lon = obj.mask.mask_lon;
                lat = obj.mask.mask_lat;
                statisticResult = nan(length(lon),length(lat));
                for i = 1:length(lon)
                    for j = 1:length(lat)
                        if obj.mask.mask(i,j) == 0
                            statisticResult(i,j) = nan;
                            continue;
                        end
                        inside_extent = [lat(j) - 1/2*obj.interval, lat(j) + 1/2*obj.interval;lon(i) - 1/2*obj.interval, lon(i) + 1/2*obj.interval]; 
                        [xbox, ybox] = EddyStatisticClass.area_box_set(inside_extent);
                        in = inpolygon(longitude, latitude, xbox, ybox);
                        statisticResult(i,j) = sum(in);
                    end
                end
            end

            function [eddyGenerationArray, eddyExtinctionArray]= GenerationAndExtinctionEddyArray(traj, type)
                eddyGenerationArray = [];
                eddyExtinctionArray = [];

                for trajI = 1:length(traj)
                    trajItem = traj(trajI);
                    eddys = EddyDatasetClass.eddyArrayAddByTraj([], trajItem, type);
                    n = length(eddys);
                    eddyGenerationArray = [eddyGenerationArray eddys(1)];
                    eddyExtinctionArray = [eddyExtinctionArray eddys(n)];
                end
            end

        end
        % % 计算涡旋属性的统计分布特征
        % function eddyDistributionstatistic(obj, varargin)
        %     idx_type = find(strcmp(varargin,'type'), 1);
        %     if ~isempty(idx_type)
        %         type = varargin{idx_type + 1};
        %         if type == -1
        %             eddyArray = EddyDatasetClass.eddyArrayAddByTraj(obj.AEEDC.traj, -1);
        %         elseif type == 1
        %             eddyArray = EddyDatasetClass.eddyArrayAddByTraj(obj.CEEDC.traj, 1);
        %         end
        %     end
        % 
        %     function distributionstatistic(eddyArray, varargin)
        %         varargin = varargin{1};
        %         idx_property = find(strcmp(varargin,'property'), 1);
        %         if ~isempty(idx_property)
        %             property = varargin{idx_property + 1};
        %         end
        %         propertyArray = [eddyArray.(property)];
        %         property_max = max(propertyArray);
        %         property_min = min(propertyArray);
        %     end
        % 
        % end
    end
    methods(Static)
        function [xbox, ybox] = area_box_set(inside_extent)
            xboundary = [inside_extent(2,:)];
            yboundary = [inside_extent(1,:)];
            xbox = xboundary([1 1 2 2 1]);
            ybox = yboundary([1 2 2 1 1]);
        end
    end
end

